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We demonstrate the suitability of tensor network techniques for describing the thermal evolution 
of lattice gauge theories. As a benchmark case, we have studied the temperature dependence of 
the chiral condensate in the Schwinger model, using matrix product operators to approximate the 
thermal equilibrium states for finite system sizes with non-zero lattice spacings. We show how these 
techniques allow for reliable extrapolations in bond dimension, step width, system size and lattice 
spacing, and for a systematic estimation and control of all error sources involved in the calculation. 

The reached values of the lattice spacing are small enough to capture the most challenging region of 
high temperatures and the final results are consistent with the analytical prediction by Sachs and 
Wipf over a broad temperature range. 
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INTRODUCTION 


Tensor network (TN) techniques have recently revealed 
their potential to study lattice gauge theories (LGT). Nu¬ 
merical studies have demonstrated that matrix product 
states (MPS) can accurately describe ground state and 
low excited levels of the Schwinger model and of 

related quantum link models 0,0, and tensor renormal¬ 
ization group methods have been used to evaluate the 
path integral 00- More generally, a framework has 
been proposed to construct gauge invariant TN states in 
higher dimensions 00. 

The TN ansatz can also be used to describe thermal 
equilibrium states. Non-zero temperature studies have 
played a major role in lattice QCD computations, estab¬ 
lishing a crossover behavior of QCD in the early universe 
[0 . Such calculations will be very important in the fu¬ 
ture for understanding QCD matter, see 14]. Different 
from most LGT calculations, tensor network methods 
work in the Hamiltonian formalism. In this approach, 
it is in principle possible to follow the complete ther¬ 
mal evolution of the system. Thus, employing Hamilto¬ 
nian calculations could allow us to map out the tempera¬ 
ture dependence of many physical quantities over a much 
broader temperature regime and in a more precise way 
than conventional Markov Chain Monte Carlo methods. 

In this work, we show that the matrix product oper¬ 
ator (MPO) ansatz can accurately describe the thermal 
equilibrium states of the Schwinger model. To this end, 
we investigate the temperature dependence of the chiral 
condensate in the continuum limit for the massless case, 


for which analytical results were provided in |15|. Our 
results are consistent with the analytical prediction at 
all (dimensionless) inverse temperatures <7/3 £ [0,6]. Our 


lattice calculation using MPO requires a series of consec¬ 
utive extrapolations. We describe how to carry out these 
steps and demonstrate that all systematic errors inher¬ 
ent to the method can be controlled and systematically 
improved. Thus, the procedure yields reliable continuum 
values and is applicable also when the exact value is com¬ 
pletely unknown. 

We use Gauss’ law to integrate out the gauge degrees 
of freedom and apply TN states to describe the fermionic 
degrees of freedom in the exact physical subspace, as in 
[3, 0. Here we demonstrate that this approach, initially 
presented in 0 , is also suitable for thermal states. Al¬ 
ternatively, one could include both fermionic and bosonic 
degrees of freedom and impose gauge symmetry on the 
tensors, as in [eim , but in that case the gauge degrees 
of freedom need to be truncated, which introduces an ad¬ 
ditional extrapolation in the procedure, which also has to 
be taken into account in the systematic errors. 


THE MODEL AND THE HAMILTONIAN SETUP 


The Schwinger model [17] or QED in 1 + 1 dimen¬ 
sions, is frequently used as a testbench for lattice cal¬ 
culations. In order to apply TN methods, we work in 
the Hamiltonian formalism (see e.g. 0 ). which implies 
we have to impose Gauss’ law as a constraint on phys¬ 
ical states. The Kogut-Susskind Hamiltonian [18| can 
be mapped to a spin system by a Jordan-Wigner trans¬ 
formation 0 . Using Gauss’ law, the gauge degrees of 
freedom can be eliminated [20], since the electric flux on 
a link is completely determined by the spin content and 
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the background field, so that the model can be written: 
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where £ is the boundary electric field (on the leftmost 
link), which can describe the background field, and the 
parameters of the model are x = , p, = , in terms 

of the lattice spacing, a, the fermion mass, m, and the 
coupling, g. 

It is then possible to use a basis \£)\ioi\. -. in-i) to 
describe the physical space. For finite systems, the value 
of i is conserved. In the following, we will consider the 
case £ = 0 and omit it from the basis. We focus on 
the temperature dependence of the chiral condensate, 
E = ('F’l ’)/g, which is the order parameter of the chiral 
symmetry breaking, and which in terms of spin operators 

reads E = ^E n (-l) ni ^- 


THE ANSATZ 

An MPS for a system of N sites with internal dimen¬ 
sion d and individual basis {|i)]q =1 is a state of the form 
W = Ei 0 ,...i w _ 1= i tr W where 

each A l k is a D-dimensional matrix, and the bond dimen¬ 
sion, 1). determines the number of free parameters in the 
ansatz [21 h 23I | . MPS are known to provide good approx¬ 
imations to ground states of local Hamiltonians in the 
gapped phase [24] and have also been successfully used 
for more general situations m- The analogous ansatz 
in the space of operators 2fJ 28 j is called matrix product 
operators (MPO), and can be used to efficiently approx- 
imate thermal states of local Hamiltonians [29|, [30(. 

To find an MPO approximation to the Gibbs state, 
p oc e~P H , a Suzuki-Trotter decomposition is applied to 
the exponential, with the Hamiltonian split into several 
terms whose exponentials are easily written or approx¬ 
imated as MPO. In the case of Hamiltonian m, it is 
convenient to split H = H e + H 0 + H z , where -H e (o) con¬ 
tains the cr+ cr~ +1 + h.c. terms for even (odd) n, and H z 
contains the mass terms and the long range cr^cr^ inter¬ 
actions. The exponentials of D e ( 0 ) can be easily written 
as exact MPO [28j. For H z , instead, the exponential can 
only be approximated. Adopting a 2nd-order Trotter ex¬ 
pansion and, for H z , a lst-order Taylor expansion, we 
can write: 
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where 5 = /3/M is the step width, and the final error for 
fixed P will be 0(6), dominated by the Taylor expansion. 


Starting from the identity operator, which corresponds 
to infinite temperature, p(P = 0), we apply successive 
Euclidean evolution steps. After each of them a trunca¬ 
tion is carried out to find an MPO approximation to the 
result. To this end, using a Choi isomorphism, the MPO 
is mapped to an MPS with local physical dimension d 2 , 
and an alternating least squares procedure is applied to 
minimize the Euclidean distance between the vectorized 
MPO for the new and evolved states. Since the trunca¬ 
tion does not preserve the positivity of the whole state, 
it is more convenient to compute p(p/2)^ p(/3/2), where 
the Trotter expansion explained above is used for each 
factor. f35l ] 


RESULTS 

To compute the chiral condensate in the infinite vol¬ 
ume and continuum limits, we approximate the ther¬ 
mal state at each g/3 £ [0, 6] over a range of values of 
x £ [4,65]. For each value, we consider various system 
sizes, with N/y/x £ [16, 24] to ensure consistent physical 
volumes, and for each of them, different step widths, 6, 
and bond dimensions. We thus need to control effects of 
successively extrapolating in D , <5, N and x. 

The limited bond dimension, D, used for each fixed 
set of values (x, N, 6) induces a systematic truncation 
error. The MPS family being complete, the results con¬ 
verge to the exact value for the given problem in the 
limit of very large D (of the order of the dimension of 
the operator space) [23|. From our finite D results, we 
estimate the final value of E as the one obtained for the 
largest achieved bond dimension, D max , and the error as 
the difference between this value and the one obtained 
from D max — 20, as illustrated in the left panel of Fig. |T] 

A second source of systematic error is the finite step 
width, 6. Although we use a second order Suzuki-Trotter 
expansion, the Taylor approximation in (J3J induces linear 
corrections, 0(5). We can thus extrapolate linearly to 
obtain the value as 5 —> 0, as illustrated by the right panel 
in Fig. [T] for selected examples. Furthermore, the Taylor 
approximation requires that the value of 5 employed in 
the calculation is below a certain threshold, to ensure 
convergence of the expansion. We find that values 6 = 
10 -6 -lCU 3 are sufficiently small for the considered values 
of x and N. 

The previous steps yield a result for each pair (x, N). 
As in 0 , 0 ], we then find the thermodynamic limit by 
fitting the results to a linear function in 1/N. The left 
panel of Fig. [2] shows how accurately this extrapolation 
fits our results for the considered values of x. 

From the infinite volume results for each lattice spac¬ 
ing at fixed g[3, we can perform the continuum extrapo¬ 
lation. As we showed in [4j , the condensate exhibits loga¬ 
rithmic corrections 0(a log(a)). Hence we try two fitting 
functions, which additionally include linear or linear and 
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quadratic corrections in a, 

fl{x) = Scant + log(x) + -^=, 

y/X y/X 

f 2 {x) = S con t + “F- log(x) + -^L + —. (3) 

\/X -y/X X 

To include the uncertainty from the choice of the func¬ 
tional form of the fitting ansatz, we finally take as central 
value the result from the fit to f\ , and as systematic error 
the difference between both. Fig. [3] demonstrates these 
fits for two very different temperatures, g/3 = 0.4 and 4.0, 
and shows that quadratic corrections will only be signif¬ 
icant at lower temperatures. Higher order corrections do 
not provide any significant improvement for the fitting 
range x G [4,65]. 

After performing for each value of g[3 all the steps de¬ 
scribed above, we obtain for the chiral condensate the 
temperature dependence shown in Fig. [5] Comparing to 
the analytical result in 0 , we find excellent agreement 
for all g/3 > 0.5. Although the central values lie very 
close to the exact results, the errors shown in Fig.0seem 
relatively large because they include the propagated er¬ 
rors from the extrapolations in D , 5 , N and x, as well as 
the systematic error from the form of the fitting ansatz 
for the continuum extrapolation. Our approach makes it 
possible to fully control all sources of uncertainties. This 
rigorous account of errors is crucial to ensure that the 
technique can be used in a general situation, for which 
no analytical results are available. 

Different kinds of errors contribute distinctly at dif¬ 
ferent temperatures. For small g/3, cut-off effects are 
enhanced, and systematic errors from the choice of the 
fitting ansatz can be an order of magnitude larger than 
other errors. Lowering the temperature, the effect be¬ 
comes smaller, while other errors grow, in particular the 
propagated error from the D —> oo extrapolation. Inter¬ 
estingly, at intermediate values, g/3 ~ 2, the slope of the 
N ^ oo extrapolations changes sign (see Fig. 0. More¬ 
over, for this region of temperatures, the continuum limit 
extrapolations from /i and f 2 are very close. For larger 
gf3, all errors grow, but the increase is faster for D and <5 
extrapolations so that at the end of our g/3 range they are 
an order of magnitude larger than the systematic error 
from the choice of the fitting ansatz. 

The approach is systematically improvable, because it 
is clear how to reduce each uncertainty. For the present 
analysis, the relatively most important error comes from 
the extrapolation in D. We chose a rather conservative 
criterion to estimate this error and it is natural to ex¬ 
pect that much more accurate results can be obtained by 
checking the convergence in bond dimension with larger 
values of D. The cost of the computation scales with 
D 3 , which makes the scan over a rather broad range of 
D feasible. It is nevertheless remarkable that the very 
accurate results presented here were obtained with a rel¬ 


atively small D < 100. This shows how adequate the 
MPO ansatz is for thermal states in this gauge theory. 

For the region of very small gf3, we find a significant de¬ 
viation from the analytical results (see Fig. 0, although 
the individual points at finite x are accurate enough, be¬ 
cause much smaller lattice spacings are required in order 
to correctly capture the asymptotic behavior. Using the 
procedure described above, it is possible to reach larger 
values of x, by incurring a higher computational cost, 
since the required system size grows as yfx to maintain 
a consistent physical volume, and correspondingly the 
threshold S that ensures convergence decreases. 

Using an alternative approximation for the exponen¬ 
tial of H, that avoids the Taylor expansion, it is however 
possible to explore the region of larger x at a lower com¬ 
putational cost. Indeed, e~ SHz can be written exactly as 
an MPO of bond dimension O(N). For systems of several 
hundreds of sites this is unpractical, but this exponential 
can be approximated by projecting out those spin config¬ 
urations that correspond to an electric flux larger than 
a certain cutoff, L c u t, on any of the links. This results 
in an MPO with bond dimension 2L cut + 1. Notice that 
this truncation is equivalent to limiting the maximum oc¬ 
cupation number for the bosonic degrees of freedom, as 
done in GLSS for pure states. However, the latter re¬ 
sults in twice as large system sizes, since bosonic degrees 
of freedom are kept explicitly, and, in principle, in local 
dimensions that scale as (2A cut + l) 2 for the additional 
sites in the MPO. 


Although the cost of applying the MPO for the pro¬ 
jected exponential is higher than that of the Taylor ap¬ 
proximation, the step width error is now 0(S 2 ), deter¬ 
mined by the second order Trotter expansion, and there 
is no threshold value for <5, which allows us to reach the 
same g/3 with fewer steps. To explore the small g/3 re¬ 
gion, we study the range x € [4,1024]. For each value 
of x, we compute different system sizes, step widths and 
A cu t values, and for each of them bond dimensions up 
to D = 160. As described above, we successively ex¬ 
trapolate 1/D —> 0 (as before), A —0 (linear in S 2 ) and 
1/N —> 0, as illustrated by Fig. [4] (left). To account for 
the additional systematic error due to the cutoff parame¬ 
ter, we can also extrapolate in L cut . However, comparing 
the results for L cut € [5,15], we observe (see Fig. 0 that 
the effect is very small, and results for L cut > 8 are com¬ 
patible within our numerical precision (inset of Fig. 0 
right), so here we present the results for L cut = 10 and 
leave the detailed analysis of the cutoff effects to a more 
technical work [311. 

This relatively small L cut allows us to study the lattice 
effects for the smallest gf3. As shown in Fig. 0 g/3 = o.i 
requires x ~ 300 or larger to obtain an accurate contin¬ 
uum extrapolation. We observe that higher order correc¬ 
tions are present and adopt as central value the result 
of the fit / 3 (x) = S C ont + ^=log(x) + ^ ^ — 


Ey/x 

37| and we estimate the error from using different fitting 
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FIG. 1: Left: convergence of condensate value with bond 
dimension, D, for g/3 = 0.4, x = 6.25, N = 40 and 8 = 
5 • 10 -4 . Shown is the deviation with respect to the final 
value, £ = 0.01508769. Looking at various 8 (inset) we find 


FIG. 3: Condensate, £, as a function of the lattice spacing for 
g/3 = 0.4 (left) and 4.0 (right). The lines show the continuum 
extrapolation using /i (solid blue) and /2 (dash-dotted red). 



FIG. 2: Infinite volume extrapolation of the condensate values 
for several values of x at gf3 = 0.4 (left) and 4.0 (right). 


ranges. Notice that to properly deal with the uncertainty 
due to the fit, we could run a statistical analysis as was 
done in {3j, but the simple estimate used here allows us 
to appreciate the relevance of reaching large x values. 

Finally, we obtain for the temperature dependence 
of the chiral condensate the improved results shown in 
Fig. [5] (right), consistent with the analytical result at 
high temperatures. 


CONCLUSION 

Using the massless Schwinger model as a benchmark, 
we have demonstrated that finite Matrix Product Opera¬ 
tors can successfully describe thermal equilibrium states 
of lattice gauge theories. We have evaluated the ther¬ 
mal evolution of the chiral condensate in this model and 


found good agreement with the analytical result [15] from 
infinite to almost zero temperature. 

The high temperature region is the hardest one to cap¬ 
ture in the continuum limit, maybe counterintuitively, 
since for typical condensed matter models it is easier 
to describe. On the other hand, it is known that lat¬ 
tice spacing effects in the high temperature region can 
be non-negligible in conventional lattice simulations. We 
have nevertheless shown that using the MPO ansatz, it is 
also possible to obtain precise results at very small lattice 


angles), 8, 10, 12 (different shapes and colors, indistinguish¬ 
able in the plot) and 15 (orange diamonds). The inset shows 
the extrapolation in Trotter parameter, 5, for N = 170 and 
L cn t = 8 — 15. Right: continuum limit for g/3 = 0.1 (left 
pointing purple triangles) and 0.2 (green diamonds). The cor¬ 
responding exact results from pji| are indicated on the vertical 
axis. The inset shows explicitly the dependence on the cutoff 
for these g/3 values at x = 55. 


spaemgs. 

Our approach offers a systematic procedure to evalu¬ 
ate and control all systematic errors in the calculation, 
namely bond dimension of the ansatz, step width of the 
Trotter expansion for the exponential operators, finite 
volume and continuum limit. Although not strictly nec¬ 
essary, a truncation of the maximum electric flux per 
link can be introduced to enhance the numerical perfor¬ 
mance. The effect of this additional cutoff parameter is 
very small, but can equally be taken into account in the 
systematic error analysis. 

All this makes the MPS and MPO ansatze most valu¬ 
able and promising tools to evaluate also other one¬ 
dimensional Hamiltonian systems relevant for gauge held 
theories. The most interesting open question is the ex¬ 
tension of these techniques to higher dimensions. 
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FIG. 5: Temperature dependence of the chiral condensate 
from data with x < 65 and exact gauge sector (red triangles), 
compared to jl5j ] (solid line). For the lowest g/3 (right), the 
results deviate from the exact ones. Using the L c u t = 10 
truncation, we reach smaller lattice spacings and recover the 
consistency with the analytical results (blue circles). 
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Numerical method 


For completeness, we describe here the details of the 
numerical method used in the paper. 

To describe the thermal equilibrium state we use the 
MPO ansatz (HQ, 


0 = 


J2 ti : ■ ■ ■ Mfrsj™- 1 ) 


{ik,jk } 


|*o • • • Uv-i) (jo • • • jjv-il • (4) 


The thermal density operator can be written as imag¬ 
inary evolution of the identity operator [261 ]. p(/3) oc 
e~^ H = e~i H te~ 2 H . By applying a Choi isomor¬ 
phism [33|, |*)(j| —l |i) ® | j), the operators can be vector¬ 


ized. The thermal state can thus be approximated as an 
MPS, by applying the imaginary time evolution operator 
corresponding to H = H ® 1 + 1 <8> H T on the initial 
vectorized identity (see Fig. [6]). As shown in Eq. 
the Hamiltonian contains non-commutative terms, so we 
approximate the exponential operator as a sequence of 
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FIG. 6: Schematic representation of the TN operations. Left: 
applying the first exponential for the first step of imaginary 
time evolution on the identity operator (which is an MPO 
with D — 1). Center: vectorizing this first step allows us to 
work with the MPS formalism. Right: after each application 
of the exponential, the result is approximated by an MPS 
with bounded bond dimension. 


MPOs, to be applied on the MPS. To this end, we use a 
2nd-order Suzuki-Trotter expansion, 


o-P H 


e -f H - e -^H Ze -6H 0e -^H Ze -^H e 


M 


(5) 


Each evolution step is thus approximated as the succes¬ 
sive action of five terms. The exponentials of the hopping 
terms, H e(o) = even (odd) 1 <Vfi + h.c.), have a sim¬ 

ple exact MPO expression with maximal bond dimension 
4 [2^, constructed as simply the product of the individ¬ 
ual exponentials of the mutually commuting two-body 
terms. The remaining term, 


JV-l 




n—0 
N—2 


4 E 0 + 0 + 2 ^^ 

n =0 L k =0 Kk 

N—2 / n \ 

EZ ( 1 + 2 a k )» 

n —0 \ k —0 / 

(even) 


-Z _z 
a k a l 
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contains long-range terms cr^o^, and, although all terms 
commute with each other, the product of individual ex¬ 
ponentials would yield a bond dimension exponentially 
large with the system size, N. 

A more efficient expression for the exponential exists 
with bond dimension that only scales linearly in N. We 
can indeed write H z as a sum of mutually commuting 
local terms, H z = J2 n h n , where, for n < N — 1, 

h n = ln + (-iro + Li (?) 


and hM-i = % [l + (— l) N ~ 1 a^_ 1 ] , L n being the electric 
flux on each link. The exponential of H z is diagonal in 
the z spin basis, and its value on a basis vector can be 
written as a product of the exponentials of each of these 
terms for the corresponding state. Since the value of L n , 
by virtue of Gauss’ law, is completely determined by the 
spin content on sites k <= n, the factor corresponding 
to a given link can be determined from the total magne¬ 
tization, ]Cfc<n°fc’ to the left of the corresponding site. 
Such information can be encoded in the virtual index of 
an MPO, which, in a chain of length N, could in prin¬ 
ciple assume values L n e [— N/2, N/2], The exponen¬ 
tial can thus be written exactly as an MPO determined 
by tensors M n whose only non-vanishing elements are 
{M") Ln _ lLn = e ~ Shn for L„ = L„_ 1 + i[(-l)” + «) ii ]. 

Such exact expression produces an MPO with bond di¬ 
mension 0{N ), which is not practical for the long chains 
involved in our study. Thus, it is convenient to truncate 
the MPO by allowing L n , i.e. the virtual index, to as¬ 
sume only bounded values \L n \ < L cut , so that the max¬ 
imum bond dimension is 2 L cut + 1. This corresponds to 
a truncation of the physical space to only those spin con¬ 
figurations for which all links have small enough electric 
flux, since the rest will be projected out when multiplying 
by the truncated exponential. 

Another, more economical approximation to the ex¬ 
ponential of H z can be achieved by a lst-order Taylor 
expansion, which can be written as an MPO with bond 
dimension 3. In this approach, the whole physical space 
is kept, so that no extrapolation in the L cut parameter is 
required. In our calculation, we use both approaches. 

The exponentials in © involve H a = H a (§) 1 +1 ® H J 
for each a = e, o, z. But since both terms in each H a 
commute, the corresponding exponential is just the ten¬ 
sor product of two exponentials, which can then be ap¬ 
plied sequentially or simultaneously. After every factor 
in ([5| is written or approximated as an MPO, the effect 
of one evolution step on a certain intermediate state, vec¬ 
torized as an MPS, can be approximated as a new MPS 
with the desired bond dimension. This is achieved by a 
global optimization (see e.g. [34] for algorithmic details), 
in which the Euclidean distance e = Hip') — 0|p)|| 2 be¬ 
tween the new MPS, | p'), and the result of the operator 
O on the original one, \p), is minimized by successively 
varying each one of the tensors, and sweeping over the 
chain until convergence (Fig. [G]). 



















